Multifractality and critical fluctuations at the Anderson transition 
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Critical fluctuations of wave functions and energy levels at the Anderson transition are studied 
for the family of the critical power-law random banded matrix ensembles. It is shown that the 
distribution functions of the inverse participation ratios (IPR) P q are scale-invariant at the critical 
point, with a power-law asymptotic tail. The IPR distribution, the multifraetal spectrum and the 
level statistics are calculated analytically in the limits of weak and strong couplings, as well as 
numerically in the full range of couplings. 



I. INTRODUCTION 

As is well known, in d > 2 dimensions a disordered 
electronic system undergoes, with increasing strength of 
disorder, a transition from the phase of extended states to 
that of localized states (Anderson transition). Another 
important realization of the Anderson critical point is 
the Quantum Hall plateau transition in a 2D system in 
strong magnetic field. One of the hallmarks of the metal- 
insulator transition is represented by strong fluctuations 
of eigenfunctions. These fluctuations can be character- 
ized by a set of inverse participation ratios (IPR) 



d d r \ip(r)\ 



2q 



(1) 



The field theory of the Anderson transition is the matrix 
non-linear er-model, in the replica [0 or supersymmet- 
ric |Q formulation. In 2 + e-dimensions with e < 1 the 
transition takes place in the weak-coupling regime, allow- 
ing for a systematic renormalization-group (RG) treat- 
ment, which yields the critical indices in the form of the 
e-expansion. In particular, Wegner H found in this way 
that the IPR show at criticality an anomalous scaling 
with respect to the system size L, 



P q cx L- T{q) 



-(g) = D q (q - 1). 



(2) 



Equation (|^) should be contrasted with the behavior of 
the IPR in a good metal (where eigenfunctions are er- 
godic), P q cx L - ^ 9-1 ), and, on the other hand, in the 
insulator (localized eigenfunctions), P q cx L°. 

The scaling (||) characterized by an infinite set of crit- 
ical exponents D q implies that the critical eigenfunction 
represents a multifraetal distribution Q] . The notion of a 
multifraetal structure was first introduced by Mandelbrot 
|^| and was later found relevant in a variety of physical 
contexts, such as the energy dissipating set in turbulence, 
strange attractors in chaotic dynamical systems, and the 
growth probability distribution in diffusion-limited ag- 
gregation; see || for a review. More recently, consid- 
erable research interest was attracted by the problem 
of Dirac fermions in a random vector potential, which 



allows for an exact calculation of the multifraetal spec- 
trum |7] |^]. The multifraetal exponents play an impor- 
tant role in recent attempts of identification of the con- 
formal theory describing the Quantum Hall plateau tran- 
sition 0,0. 

During the last decade, multifractality of critical eigen- 
functions at the Anderson transition has been a subject 
of intensive numerical studies, see Refs. @Jl| and ref- 
erences therein. Among all the multifraetal dimensions, 
D2 plays the most prominent role, since it determines 
the spatial dispersion of the diffusion coefficient at the 
mobility edge [14j|. 

It should be stressed, however, that Wegner's result (|^) 
refers to an ensemble-averaged IPR. On the other hand, 
it is well known that in disordered systems mesoscopic 
fluctuations from one realization of disorder to another 
may be very strong. As a result, an average value of some 
quantity may not provide sufficient information, and one 
has to speak about the corresponding distribution func- 
tion. This poses the question of the statistics of the IPR's 
P q at criticality, which is a central issue of the present 
article. 

Let us first remind the reader of the existent analyt- 
ical results concerning the IPR fluctuations. While the 
direct analytical study of the Anderson transition in 3D 
is not feasible because of the lack of a small parameter, 
statistics of energy levels and eigenfunctions in a metallic 
mesoscopic sample (dimensionless conductance g » 1) 
can be studied systematically in the framework of the 
supersymmetry method; see |l5[ for a review. Within 
this approach, the IPR fluctuations were studied recently 
p6| , p7 15 1 . In particular, the 2D geometry was consid- 
ered, which, while not being a true Anderson transition 
point, shows many features of criticality, in view of the 
exponentially large value of the localization length. It 
was found that the distribution function of the IPR P q 
normalized to its average value (P q ) has a scale invariant 
form. In particular, the relative variance of this distribu- 
tion (characterizing its relative width) reads 



var(P g ) _ Cq 2 (q-l) 



(P q ) 2 



P 2 g 2 



(3) 



where C ~ 1 is a numerical coefficient determined by the 



1 



sample shape (and the boundary conditions), and [3—1 
or 2 for the case of unbroken (resp. broken) time re- 
versal symmetry. It is assumed here that the index q is 
not too large, q 2 <C fiirg. These findings motivated the 
conjecture that the IPR distribution at criticality 
has in general a universal form, i.e. that the distribu- 
tion function V (P q / Pq VP ) is independent of the size L in 
the limit L — > oo. Here P q yp is a typical value of the 
IPR, which can be defined e.g. as a median |L8| of the 
distribution V(P q ). Normalization of P q by its average 
value (P q ) (rather than by the typical value Pq YP ) would 
restrict generality of the statement; see the discussion 
below. Practically speaking, the conjecture of Ref. |l^] 
is that the distribution function of the IPR logarithm, 
VQxiPq) simply shifts along the x-axis with changing L. 

The applicability of these results to a generic Anderson 
transition point has been questioned recently in Ref. [|l9| . 
Indeed, the 2D metal represents only an "almost crit- 
ical" point, and the consideration is restricted to the 
weak disorder limit j> 1 (weak coupling regime in the 
field-theoretical language), while all the realistic metal- 
insulator transitions (conventional Anderson transition 
in 3D, quantum Hall transition etc.) take place in the 
regime of strong coupling. It was proposed in (on the 
basis of numerical simulations for the 3D tight-binding 
model) that the fractal dimension D2 is not a well de- 
fined quantity, but rather shows universal fluctuations 
characterized by some distribution function V(D^) of a 
width of order unity. 

To explore the IPR fluctuations (and also the level 
statistics, see below) at criticality in the full range from 
weak to strong coupling, we consider the power-law ran- 
dom banded matrix (PRBM) ensemble. The model is 
defined 20 1 as the ensemble of random Hermitean N x N 
matrices H (real for (3 = 1 or complex for (3 = 2). The 
matrix elements Hij are independently distributed Gaus- 
sian variables with zero mean (Ha) = and the variance 



(|^| 2 )=a 2 (|i-j|) 



where a(r) is given by 
a 2 (r) = 



1 + (r/6) 2 



(4) 



(5) 



At a = 1 the model undergoes an Anderson transition 
from the localized (a > 1) to the delocalized (a < 1) 
phase. We concentrate below on the critical value a = 1, 
when a(r) falls down as a(r) oc 1/r at r 3> b. Morejare- 
cisely, we will study the periodic generalization of 1 



a 2 (r) 



1 sin 2 (7rr/A0 
fo 2 (tt/A) 2 



(6) 



(an analog of the periodic boundary conditions), which 
allows us to get rid of the boundary effects. 

In a straightforward interpretation, the PRBM model 
describes a ID sample with random long-range hopping, 



the hopping amplitude decaying as l/r a with the length 
of the hop. Also, such an ensemble arises as an effec- 
tive description in a number of physical contexts, such 
as the quantum Fermi accelerator [[2l) , the derealization 
of impurity-induced quasiparticle states in a 2D <i-wave 
superconductor [p2[ , the scattering by a Coulomb center 
in an integrable billiard |22[ |, the motion of two inter- 
acting particles in a ID random potential [pi) , and the 
quantum chaos in a billiard with a non-analytic bound- 
ary pq] . Very recently, a connection between the level 
statistics of the PRBM model at b 3> 1 and the corre- 
lations in the Luttingcr liquid at finite temperature has 
been established |26| ]. 

At a = 1 the PRBM model was found to be criti- 
cal for arbitrary value of 6; it shows all the key features 
of the Anderson critical point, including multifractality 
of eigenfunctions and non-trivial spectral compressibil- 
ity p0 15 1 . The important property of the ensemble is 
the existence of the parameter b which labels the criti- 
cal point: Eqs. (|4|), (^3j) define a whole family of critical 
theories parametrized by b [^tJ. This is in full analogy 
with the family of the conventional Anderson transition 
critical points parametrized by the spatial dimensionality 
2 < d < 00. The limit b 3> 1 is analogous to d = 2 + e 
with e <C 1; it allows for a systematic analytical treat- 
ment via the mapping onto a supermatrix a-model and 
the weak-coupling expansion |2^,^]. The opposite limit 
!)< 1 corresponds to d ^S> 1, where the transition takes 
place in the strong disorder (strong coupling) regime. As 
we demonstrate below, it is also accessible to an ana- 
lytical treatment using the RG method introduced by 
Levitov |2£| |. Let us also note a similarity with confor- 
mal models proposed recently as candidate theo- 

ries of the Quantum Hall critical points, which are also 
parametrized by a continuously changing marginal cou- 
pling constant. 

In this paper we will combine the analytical study of 
the eigenfunction and energy level statistics in the two 
limits b 1 and b <C 1 with numerical simulations in 
the full range of b. The feasibility of the systematic an- 
alytical treatment of both regimes, weak-coupling and 
strong-coupling, as well as of the numerical simulations 
in a broad range of couplings, makes the PRBM ensem- 
ble a unique laboratory for studying general features of 
the Anderson transition. 

As has been already mentioned, we will study not only 
the eigenfunction fluctuations, but also the energy level 
statistics. It has been understood |]30|-^] that a scale- 
invariant level statistics is an intrinsic feature of the An- 
derson critical point. In particular, the critical level cor- 
relations are characterized by a non-trivial value of the 
spectral compressibility < x < 1 (intermediate be- 
tween y = in a metal and x = 1 m an insulator) 
pfj| , ^3 34 1 . While the level correlation function itself is 
shape-dependent, the value of x is a fully universal at- 
tribute of the critical theory (like critical indices). 

Supplementing again the analytical study at b ^> 1 and 
b <C 1 by numerical simulations, we are able to calculate 



2 



the two-level correlation function i?2(s) and the spectral 
compressibility \ m the whole range of b. Our interest 
to the critical spectral statistics was additionally moti- 
vated by a recent paper J3q] , where it was argued that 
the following exact relation between \ an d P>i holds: 



X 



d-D 2 
2d 



(7) 



According to (0) , the spectral compressibility should tend 
to 1/2 in the limit D2 — > (very sparse multifractal), and 
not to the Poisson value x = 1. The derivation of (Q) is 
based, however, on a certain assumption of the decou- 
pling of the energy level and eigenfunction correlations, 
which is difficult to verify directly. It is important, there- 
fore, to check the validity of the result (0), and such an 
opportunity is provided by the PRBM model. 

The structure of the article is as follows. Sect, g is 
devoted to the IPR statistics and the multifractal spec- 
trum of the PRBM model. In Sect. Ill we study the 



two-level correlation function and the spectral compress- 
ibility. Sect. [V summarizes our findings. Some of the 



results of this work have been presented in a brief form 
in the Letter [06| . 



II. MULTIFRACTAL SPECTRUM AND THE IPR 
STATISTICS 

We find it convenient to organize this section in the 
following way. 



We begin by formulating in Sect. |II A 
our main results conc erning the IPR statistics. Then we 
present, in Sect 



II B and Sect. II 



the analytical cal- 
culations in the limits of b <C 1 and i) > 1, respectively. 
The numerical simulations (which have been performed 
for (3 — 1) fully support the analytical findings, ascertain- 
ing that the approximations made in the RG treatment 
are justified. Also, they allow us to explore the interme- 
diate range of b ~ 1. A discussion of finite-size effects 
playing an important role in the analy sis of the scaling 
of the IPR distribution, is given in Sec. II D. 



A. General considerations and a summary of the 
results 

For further needs we define two sets of fractal expo- 
nents characterizing the scaling of the average IPR (P q ) 
and of the typical IPR P q yp , respectively: 



(P q ) oc N- f («) , f(q) 
P*yP oc N- T ^ , r(q) 



D q (q-l) ; 
D q (q-1) . 



(8) 
(9) 



Note that we consider q > only; negative q are outside 
the range of applicability of our analytical methods. As 
has been already mentioned, P q yp can be defined as a 
median of the distribution V(P q ) [JL8| ; an alternative def- 
inition can be P q yp — exp(lnPg). Obviously, an informa- 
tion about the IPR distribution function V(P q ) is needed 



in order to judge whether the exponents f(q) and r(q) 
are equal to each other or not. As we will demonstrate 
below, in the limit of large system size N the distribu- 
tion V(P q /P q yp ) becomes independent of N. An impor- 
tant property of this scale-invariant IPR distribution is 
its power-law "tail" at large P q , 



v{p q /Pl yp ) cx {Pjp^y 



p q » P* yp . (10) 



Of course, the far tail of this distribution becomes in- 
creasingly better developed with increasing N . In other 
words, the point where the distribution deviates from its 
limiting scale-invariant form moves to infinity as N in- 
creases. 

It is clear that the relation between r(q) and f(q) de- 
pends crucially on whether the power-law exponent x q is 
smaller or larger than unity. If x q > 1, the two definitions 
of the fractal exponents are identical, r(q) = f(q). This 
situation will be shown to occur at not too large values 
of q; in particular, xi > 1 at any b. However, with in- 
creasing q the value of x q decreases. Once it drops below 
unity, the average (P q ) starts to be determined by the 
upper cut-off of the power-law "tail" , which depends on 
the system size. As a result, (P q ) shows scaling with an 
exponent f q different from r q . In this situation the aver- 
age value (P q ) is not representative and is determined by 
rare realizations of disorder. 

The connection between f(q) and r(q) in the regime 
x q < 1 can be elucidated best via introducing the singu- 
larity spectrum f(a), which is the conventional way of 
analyzing multifractal distributions ||. To this end, let 
us note that the average IPR's (P q ) are (up to a multipli- 
cation by N) the moments of the distribution 'PdV' 2 !) of 
the eigenfunction intensities. The behavior (j|) of the mo- 
ments corresponds to the intensity distribution function 
of the form 



W 2 I) 



1 

W 



(11) 



Indeed, calculating the moments (l^ 29 !) with the distri- 
bution function (11), one finds 



(P q ) = N(\^\) 



daN- qa+ K a) 



(12) 



where we have introduced a = — ln|^/> \/hxN. Evalu- 
ation of the integral by the saddle-point method repro- 
duces the result flq), with the exponent f (q) related to the 
singularity spectrum f(a) via the Legendre transform 



f(q) =qa- f(a) 



q = f'(a) 



(13) 



It is not difficult to see that the condition x q = 1 is 
equivalent to f(a) = 0. Indeed, both conditions x q < 1 
and f(a) < characterize the situation when the average 
value (P q ) is not representative and is determined by rare 
realizations of disorder. On a more formal level, this can 
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be derived from the formula relating x „ an d the fractal 
exponents, see Eq. ( |l6| ) below and Sec. |ll C| . 

We further denote the value of a determined by f(a) = 
as a_, and the corresponding value of q as q c (clearly, 



both 



and q c depend on b). The value of r q in the 



region q > q c can be found by observing that P q yp can 
be written in the form similar to ([L2l 



ptyp 

9 



(14) 



/(q)>0 



The restriction on the integration range removes from 
consideration the rare events of such large values of \ifi\ 2 
which can be found only in a small fraction (~ N^ a ^ 
with f(a) < 0) of all eigenfunctions. Since for q > q c the 
saddle-point a < a_ is outside the integration domain, 
the integral (jl4|) is determined in this case by the bound- 
ary a_ of the integration range, yielding (see a related 
discussion in H) 



r(q) = qa. 



q> q c 



(15) 



The value of the power-law-tail index x q is related to 
the fractal exponents as follows: 



x q r{q) = f(qx g ) 



(16) 



To be precise, we were able to derive Eq. ( |l6| ) for all 
q in the limit b <C 1, as well as for integer values of 
x q = 1,2,... at arbitrary b. We expect, however, that 
this relation is generally valid. 

According to what has been said above, the curve q c (b) 
separates the regions with the two different types of the 
multifractal behavior: at q < q c {b) we have x q > 1 and 
r(q) = f(g), while at q > q c (b) the tail index x q < 1 and 
r(q) is different from f(q) and given by (|l|). We have 
calculated the asymptotic form of the "phase boundary" 
q c (b) in both limits b ^> 1 and b <C 1, 



q c (b) 



(27T/36) 1 / 2 
2.4056 , 



6> 1 
6< 1 



(17) 



Notice that q = 2 always belongs to the low-g phase, 
i.e. t (2) = f(2) for all b. For q > q c (b) we find from 
Eqs. ([if]), ( pif ) that x q = q c (b)/q, while in the opposite 
regime q < q c (b) the value of x q is determined by the 
form of the function f(q). In particular, at b > 1 wc 
have 



f(q) = (q- l)D q ~ (q - 1)(1 - q/2(3irb) , 



(18) 



yielding x q ~ 2(3iTb/q 2 for g < (2/37T6) 1 / 2 . In the other 
limit, b <C 1, the function f(q) has a somewhat more 
complicated form 



r(«) 



46 r(g-l/2) 



(19) 



A r(«-i) ' 

and Eq. (|l^) does not seem to have an analytical so- 
lution for x q . However, for the particularly important 
case q — 2 we find X2 — 3/2, while all higher inte- 
ger q = 3,4, .. . are already above the phase boundary 
<7 C (6 « 1) ~ 2.4. 



B. Regime b > 1 



The regime b 3> 1 can be studied via the mapping onto 
the supermatrix tr-model ^,^5[ . The a- model action has 
in momentum space the form 



S[Q] =/3Str 



1 



dk . . i-Kvuj 
. . ^\k\QkQ-k —Qo^ 



. (20) 



where Q fc = ^ r e lfcr Q(r) and Q(r) is a 4 x 4 (/3 = 2) 
or 8 x 8 (/? = 1) supermatrix field constrained by 
Q 2 (r) = 1 and belonging to a coset space with the origin 
A = diag(l, —1). Furthermore, v is the density of states 
given by the Wigner semicircle law 



HE) 



2ir 2 b 



(4tt& - E 



2U/2 



\E\ < 2Virb 



(21) 



and t <C 1 is the coupling constant, 



7 = 7^"rh 
t 4 



2,2 



- 1 



4:7Tb 



(22) 



For a system of finite size iV with the periodic general- 
ization (^) of the 1/r decay law of a(r) the /c-integration 
in (EOh is replaced by summation in the usual way: 



dk 
2^ 



1 

N 



E 

=2Tm/N\ n=0,±l,±2,. 



The eigenfunction statistics can now be studied via the 
same methods as for conventional metallic samples. The 
main difference between the action (||) and that of the 
diffusive tr-model is in the replacement of the diffusion 
operator ^Dk 2 by \\k\. Consequently, all calculations 
within the weak coupling expansion of the u-model are 
generalized to the PRBM case by substituting II(fc) = 
f/8|fc| for the diffusion propagator II(fc) = 1/nvDk 2 . In 
particular, calculating the average IPR 
the following result for the fractal dimensions D q 



P q ), one finds 



t 47T/3 . , 

D** 1 -***' q< t- (23) 



Fig. ^] shows that this result is in good agreement with 
numerical simulations. 
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FIG. 1. Fractal dimension D2 as a function of the parame- 
ter b of the PRBM ensemble. The data points are the results 
of the numerical simulations, while the lines represent the 
b ^> 1 and i « 1 analytical asymptotics, D2 = 1 — l/nb 
[Eq. @] and D 2 = 2b [Eq. @]. 

Similar to the case of 2D |3 or (2 + e) dimensions §, 
Eq. ( p3| ) describes weak multifractality: the deviation of 
the fractal exponents from the "normal" value d = 1 is 
proportional to the small parameter t. The Legendre 
transform of f(q) = (g— 1)-D 9 yields the singularity spec- 
trum (for definiteness, we concentrate on the band center 
E = 0, where t = 4/6) 



/(«) 



(qp - a) 2 
4(a - 1) 



ao = 1 



2/3tt6 



(24) 



which crosses the x-axis at the point 



a_ 



(2/3tt6) 1 /2 



(25) 



corresponding to q c (b) = (2/37T6) 1 / 2 . 

In Fig. we confront our analytical findings with data 
from numerical simulations. At b — 4 the parabola repre- 
sents the numerical data well up to q ~ 8. The deviations 
from the asymptotic (parabolic) form are much more pro- 
nounced at b = 1. These deviations are a precursor of 
the crossover to the small-6 regime (Sec. |Xc]), where the 
parabolic approximation breaks down completely. The 
sign of the deviations (downwards) is consistent with the 
fact that at b = l/2ir the parabolic approximation would 
predict a_ = 0, while we expect a_ > for all b, in view 
of the absence of localized states. 




-0.8 



FIG. 2. Multifractal spectrum /(a) for b = 1 (□) and 6 = 4 
(o). Solid line indicates the parabolic approximation Eq. 
Inset: exponent f(q) (□) and r(q) (■) for 6=1. 

Now we turn to the IPR fluctuations, which are found 
by generalizing the results obtained for metallic samples 
]l6| , p7 Jl5 1 . In particular, the IPR variance is given for 
q « q c (b) by @ 



var (A) _ ±2, _ J_ V- U 2 (k) 



1 q 2 {q-lf 



24/3 2 



b 2 



(26) 



where the fc-summation goes over the non-zero harmon- 
ics k = 2-Kj/N with j = ±1, ±2, . . .. Equation ( p6[ ) is the 
PRBM counterpart of formula (||) for 2D metallic sys- 
tems. The higher moments of the IPR distribution were 
studied by Prigodin and Altshuler [|l7| (see also Jul); 
generalizing these results, we find that the irreducible 
moments (cumulants) of the order 2 < n <C 2ir(3b/q 2 are 
given by 



<U?» _ (n-l)! 



9(9-1) 



— V IL n (k) 



"-'"I'^'l 



(27) 



where C\n) is the Riemann (^-function. Defining in anal- 
ogy with |l7j 



P = 



(Pa 



- 1 



2tt/36 



(28) 



we have for the cumulants of P 
(P)=0; 

((P n )) = (n-l)\an) = K n , 



n = 2,3, . 



(29) 
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This allows us to restore the corresponding distribution 
function: 



nr(q) = f(nq), 



(33) 



V(P) = 



d s isP 

7T e ex P 

^i e i S (p+c )r(1+ . s) 
2n 



n! 

.n=2 



-P-C 



exp(— e 



-P-Cs 



(30) 



where C ~ 0.5772 is the Euler constant. The restriction 
on n given above (|2^) implies that Eq. (|3^) is valid for 
P q /{P q ) - 1 « 1. 

The similarity with the 2D metallic regime extends also 
to the asymptotic behavior of the distribution. Specif- 
ically, at P q j (P q ) — 1 ~ 1 the exponential falloff j30| ) 
crosses over to a power-law tail (see [[l5| for the discus- 
sion of this tail in 2D) 

P{P q ) ~ (PJiP,))- 1 -** ■ (31) 
To calculate x q , we consider the moments 

(p?)= E mn)\ 2q ...mrn)\ 2q . (32) 
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FIG. 3. Distribution function V{P q ) at q = 2 (o), 4 (□), 
and 6 (o) at 6 = 4 for systems of size N = 4096. The solid line 
represents the analytical result Eq. (pic|). The scattering of the 
data at small values of V is due to statistical noise. (Number 
of matrices in the ensemble: 428) Inset: Asymptotic of V(Pi). 
Dashed line indicates power law with exponent xa = 1.7. 

For n not too large the sum is dominated by the main 
part of the total configuration space, with all points n ly- 
ing far from each other, so that (P™) - (P q ) n ~ N~ nT ^\ 
In contrast, for sufficiently large n the integral is dom- 
inated by the contribution from r\ « r-i « . . . rs r n , 
yielding (PI 1 ) ~ N- f ^\ Therefore, if 



we have the marginal situation, which implies that 
V(Pq) oc P^ 1- ™. This completes the derivation of 
Eq. ( |l6| ) in the range q < q c (b). Using now ( |T6| ) in com- 
bination with Eq. (Eq), we find 



2ndb , 
—5- , q 2 < 27T/36 
q 2 



(34) 



Note that an analogous consideration in the 2D case 
yields x q — 2irf3g/q 2 , in full agreement with the result 
of the optimum-fluctuation method |l5[ ] . 

We have therefore presented an explicit calculation of 
the IPR distribution function at b 1 and q < q c (b). 
The distribution function is scale-invariant and has the 
form (p0|), (||) at P q /(P q ) - K 1 and the power-law tail 
©, (|f) at Pg/(P g ) - 1 > 1. Fig. | shows results of the 
numerical simulations for the distribution of the IPR's 
P q with q = 2, 4, and 6 at & = 4 (the corresponding value 
of q c being g c = (871-) 1 / 2 ~ 5). It is seen that at q = 2 
the analytical formula (^0|) nicely describes the "main 
body" of the distribution, with the upward deviations 
at large P indicating the crossover to the power-law tail 
(pl|). The asymptotic behavior ( pl| ) is outside the reach 
of our numerical simulations for q — 2, however, since the 
condition of its validity P ^> 2ir(3/q(q — 1) ~ 12.5 corre- 
sponds to very small values of the distribution function 
P(P) -C 10 -5 , and its clear resolution would require a 
much larger statistical ensemble. The situation changes, 
however, with increasing q (see the data for q — 4, 6 in 
Fig. ||). Equation ( |30| ) gets inapplicable (since the con- 
dition of its validity q <C q c is not met anymore), and 
the power-law asymptotic behavior (pfi) becomes clearly 
seen. In particular, the inset of Fig. ||~shows the tail for 
q — 4; the extracted value of the index 24 ~ 1.7 is in 
good agreement with the prediction of the 6^1 theory, 

X4 = 7r/2. 

In conclusion of this subsection, we comment on the 
notion of the termination of the multifractal spectrum 
(pi|), which has been discussed in the literature on dis- 
ordered Dirac fermions Q. It is important to realize 
that, in the present context, there are two types of such 
termination, depending on whether one studies P* yp or 
(P q ). In the former case, the relevant values of a are 
those where /(a) > 0, so that the singularity spectrum 
/(a) effectively terminates at a_ (which corresponds to 
q = g c (b) ~ (27T/36) 1 / 2 ). For q > q c (b) r(q) is given by 
Eq. (|l5|), so that the fractal exponent D q = qot-/(q — 1) 
saturates as D q — > a_ in the limit q — * 00. In contrast, 
if the exponent f(q) describing the scaling of the average 
(P q ) is studied, then the behavior ( p3|) continues up to 
q sa 7T/36, which corresponds to a = 0. This type of ter- 
mination (which takes place at q parametrically much 
larger than q c (b)) has a physically transparent origin: 
since a = — In \ ^ 2 \/ In TV and in view of the wave function 
normalization, the allowed values of a are restricted by 
a > 0. More detailed discussion of the behavior of f(q) 



G 



in the vicinity of q = irfib (i.e. the precise form of this 
termination) is outside the scope of the present article. 



C. Regime b < 1 

In the case b <§; 1 the problem can be studied via the 
renormalization-group method of Levitov |2§| , |38|| . The 
idea of the method is as follows. One starts from the 
diagonal part of the matrix H, each eigenstate being lo- 
calized on a single site. Then one includes into consider- 
ation non-diagonal matrix elements Hij with d(i,j) = 1, 
where d(i,j) is the distance between the sites i and j 
with periodic boundary conditions taken into account, 



d(i,j) = mm{|£- j\, N - \i - j\} 



(35) 



Now one argues that most of these matrix elements are 
essentially irrelevant, since their typical value is ~ b, 
while the energy difference \Ei — Ej\ is typically of order 
unity. Only with a small probability (~ b) is \Ei — Ej\ 
also of the order of b, so that the matrix element mixes 
strongly the two states. Following Levitov, we will say 
that these two states are in resonance. In this case one 
is led to consider a two-level problem 



two— level — 



Ei V 
V E 2 



(36) 



The corresponding eigenfunctions and eigenenergies are 



(+) _ 



COS ( 

sin ( 



</;(-) = 



— sm( 
cos 9 



where 



E+ = 



tan# 



E x +E 2 



± \V\^1+T 2 



T — , OJ = El — En 

2V ' 



(37) 
(38) 

(39) 
(40) 



In the next RG step the matrix elements with 
d(i,j) — 2 are taken into account, then those with 
d(i,j) — 3, and so forth until d(i,j) = N/2. Each time a 
resonance is encountered, the Hamiltonian is re-expressed 
in terms of the new states. Since the probability of a res- 
onance at a distance r = d(i,j) is ~ b/r, the typical scale 
T2 at which a resonance state formed at a scale r% will be 
again in resonance satisfies 



1 r2 1 
In , 

7*1 



(41) 



so that ri is much larger than r\. Therefore, when con- 
sidering the resonant two- level system at the scale r 2 , one 
can treat the ri-resonance state as point-like. Further- 
more, it is easy to see that the Gaussian statistics of the 



matrix element coupling the states on the scale r 2 is not 
affected by the transformation to the new basis induced 
by the T\ -resonance. 

Now we consider the evolution of the IPR distribution 
with the distance r; we will denote the corresponding dis- 
tribution function as f(P q , r). When a resonance occurs, 

two states with IPR's Pq and Pq 2 ^ are replaced by two 
new states with the IPR's 



P q + ) = P q lj «'* 



2q , 



PP Sill 



2«/ . 



pi-) = p« S in 2 « 9 + PW cos 2 « 9 . 
We thus have for real matrices {(3 = 1) 



(42) 



d 



f(P q ,r)=2u duj 



dV 



-V 2 f 2 /2b 2 



x[-/(P„r) + J dPWdPWf(pW,r)f{pW,r) 
xS{P q - cos 29 9 - P ? (2) sin 29 9)} , (43) 

where v — —^=er E I" 1 is the density of states and 



N nr , „ 

r = — sin — . (44) 

The difference between f and r is irrelevant for the 
present calculation, since the r-integral will be of loga- 
rithmic nature and thus dominated by r <ti N. However, 
this difference is important for the calculation of the level 
correlation function below (Sect. III). Transforming the 
integration measure according to 



duj = 2Vdr , dr 



1 



2 sin 2 9 cos 2 



-d0 



(45) 



calculating the 1^-integral, and specializing on the band 
center (E = 0) for definiteness, we reduce the evolution 
equation (|||) to the form 







dhir 



2b I"*/ 2 

JO 



d9 



sin 9 cos 2 9 



x[-f(P g ,r) + J dpUdpWf( P W,r)f(p( 2 \r) 
x5(P q - P (1) cos 29 9 - P (2) sin 29 9)} . (46) 



Eq. (^) is a kind of kinetic equation (in the fictitious 
time t = b In r) , with the two terms in the square brackets 
describing the scattering-out and scattering-in processes, 
respectively. 

Figure |4| shows the results of the numerical integra- 
tion of Eq. ( pi| ) for q — 2 with the initial condition 
/(P2) = <5(p2 — 1) at t = Q. It is seen that at sufficiently 
large t the distribution of In P2 acquires a limiting form, 
shifting with t without changing its shape. This conclu- 
sion of scale-invariance of the IPR distribution will be 
corroborated below by analytical arguments. 
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FIG. 4. Flow of the distribution of InP? calculated from 
the kinetic equation ( |l6| ) at t — blnr = 1.2 . . . 1.7 (from right 
to left). The oscillations near lnP2 = —1.5 are numerical 
artifacts due to rounding errors. 

Turning to the theoretical analysis, we consider first 
the average value (P q ). Multiplying Eq. (flfj) by P q and 
then integrating over P 9 , we get the evolution equation 
for (P a ) 



with 



f(q) = - 



d(P q 
d In? 



7T/2 



-2bf(q)(P q ) 



(47) 



d.6 



sin 2 9 cos 2 



-(l-cos 2<? 6>-sm 



2<l , 



2 r( g -i/2) 
r(g-i) 
i r(2 g -i) 

22?-3 r(g)r(g - 1) ' 



(48) 



We assumed here that q > 1/2, which is the condition of 
the existence of the integral in Eq. ( |4g| ) . For smaller q the 
resonance approximation breaks down. Integrating ( pl^ ) 
from r — 1 to r ~ N, we find the multifractal behavior 
(P q ) ~ N~ f (^ with the exponents 

f(q) = 2bf(q) . (49) 

The function T(q) is shown in Fig. ||. Its asymptotics are 



f(q)± 
f(q) 



" 77( 9 -l/2) ' 

_2_ ff l/2 



5-1/2 
q > 1 . 



(50) 
(51) 



Legendre transformation of (^9|) produces the /(a)- 
spectrum of the form 



where F(A) is the Legendre transform of T(q). The func- 
tion F(A) is shown in the inset of Fig. ^, its asymptotics 
are 



F(A) 
F(A) ~ 4 



i 

ttA 
A 

2 ' 



A-> ; 
A — > oo . 



Furthermore, it changes sign at 
sponding to g c ~ 2.4056. 



(53) 
(54) 

0.5104, corre- 




FIG. 5. Universal function T(q) characterizing the expo- 
nents t (q) via r(q) = 2bT(q) at b <C 1. Dashed line indicates 
the pole position. Inset: Legendre transform F(A) describing 
the multifractal spectrum via f(a) — 2bF(a/2b). 

These analytical findings are fully supported by nu- 
merical simulations as can be seen from Fig. H. 




f(a) = 2bF(A) ; A = a/2b , 



(52) 



FIG. 6. Multifractal spectrum /(a) for b = 0.25 (O) and 
b — 0.1 (A). Inset: exponent f(q). Dashed and dotted lines 
indicate the analytical results Eqs. (p2j) and (M. 
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We return now to the IPR distribution function. 
The scale invariance of the limiting distribution has 
been already demonstrated via the numerical solution of 
Eq. (^), see Fig. ||. To show this also analytically, we 
make the Ansatz 



f(P q ,r)=r^>f (P q r^>) . 
Substituting ( |55| ) into (f46|), we get the equation 



(55) 



2b r' 2 



d6 



r(q)[fo(P q ) + P q fo{P q )] 

tt J sin ft cos z 

x[-/o(P g ) + 1 dPWdpWf (pU)f (pM) 

x 5(P q - cos 2 « 9 - P g (2) sin 29 0)} . (56) 

The fact that the scale r has dropped out from Eq. J56| ) 
implies the consistency of the Ansatz ([35]) for the fixed- 
point distribution. To demonstrate the significance of 
this statement, we make a more general Ansatz for the 
limiting distribution 



f(P q ,r) = 



1 f\n(P q r T ^) 



Pqo{r) V 



a(r) 



(57) 



which allows for a change of the width of the distribution 
of hiP q with r, characterized by a function a(r) (note 
that cr(r) is defined up to a constant factor, which can 
be absorbed into the definition of the function go)- At 
a(r) — 1 this reduces to our earlier Ansatz (jig), while 
at cr(r) = lnr we get the form proposed in Ref. [fl9|| . 
Substituting ( |57| ) in (^6[), we find that the r-dependence 
drops out of the resulting equation for g$ if and only 
if cr(r) = const. This means inconsistency of the Ansatz 
© with a non-constant <r(r) and, in particular, excludes 
the possibility of a universal distribution of fractal expo- 
nents [<j(r) = lnr] advocated in Ref. [lS| . 

We turn now to the power-law tail of this scale- 
invariant distribution, fo(P q ) ~ P _:r5_1 . In order to 
calculate the index x q , we consider Eq. (|5^) in the limit 

P q » 1. It is easy to see that the integral J dP q ^dP q 2 ^ . . . 



in the r.h.s. of (|5q) is dominated by the region P c 



(i) 



(2) 
? 

(1) 



5( 2 ) 



p 



1 (or vice versa), the contribution of the region 
^ P q being suppressed by an additional fac- 
Furthermore, when Pq ~ P q and Pq 2 ^ ~ 1, 
in the argument of the <5-function. 

become then trivial, 



tor of P q 

we can neglect Pq 2 



(i) 



The integrals over P q 
and Eq. (1561) reduces to 



and P g (2) 



T (9> 9 = — / 



rfft 



sm 2 6> cos 2 I 



(1 - sin 29 *' ft - cos : 



2ga; g 



(58) 

Comparing this with (fl8|), (^), we see that the r.h.s. 
of (H|) is nothing but f(qx q ), so that Eq. j5q) can be 
rewritten in the form (|T^). 



We analyze now Eq. ( [l6] ) in the regimes q < q c and 
q > q c . In the case q < q c we expect x q > 1 and 
r(g) = f(g). The latter statement can be directly 
proven by applying the operation J dP q P q . . . to Eq. ([56]). 
[Clearly, this proof breaks down for q > q c because of 
the divergence of the integral J dP q P q fo(P q ).] Graph- 
ical interpretation of Eq. ( |l6| ) for q < q c is shown in 
Fig. [?]; its solution x q > 1 decreases with increasing q, 
reaching unity at q — q c , as expected. For the most 
frequently studied case q — 2 ("conventional" IPR) we 
find Xi = 3/2. As to the q > q c reg ime, we have then 
T (q) = Q a — ) an d the solution of ( |l6|) has a very simple 
form 



q 



q> q c 



(59) 




FIG. 7. Graphical interpretation of Eq. ( |l6| ) for q < q c . 

Let us remind the reader that up to now we considered 
in this subsection the ensemble of real matrices (ft = 1). 
However, all the above considerations are also applicable 
to the case ft = 2, with a minor modification. Specifically, 
the measure of the ^-integration in Eq. (^3|) should be 
modified: 



dV- 



-V 2 f 2 /2b 2 



2f 2 

dVV-^e 
tr 



V 2 f 2 /b 2 



(60) 



This leads, after the ^-integration, to the replacement 
of b by T^b. With this substitution, all results of this 
subsection remain valid for ft = 2. 



D. Finite-size effects in the scaling of the IPR 
distribution 

Since in reality one always has to deal with systems 
of a finite size, the understanding of finite-size effects 
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is important for an accurate interpretation of numerical 
data. In Figs. || we show the evolution of the distribution 
■p(lni-2) with N for three values of 6, representative of 
the small-6, the large-6, and the crossover regimes. 

In the case b -C 1, the evolution of the IPR distribu- 
tion with TV is governed by the "time" t — 6 In AT, so that 
t 3> 1 is the condition of the proximity to the fixed point. 
Therefore, at small b one needs exponentially large values 
of N in order to reach the limiting distribution. (Note 
that this is not true for the average (P q ), the evolution of 
which is governed by Eq. (Et|), implying a much weaker 
condition N ^> 1 for the scaling regime.) The logarithmi- 
cally slow approach to the limiting distribution is clearly 
seen in Figs. [§L ||. 
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FIG. 8. Evolution of the distribution V(\nP 2 ) for a) b = 1, 
b) 6 = 4 and c) 6 = 0.1 with the system size TV (from left to 
right: N = 4096, 2048, 1024, 512, (256)) 
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FIG. 9. Maximum value of VQnPz) for b — 0.1 as a func- 
tion of the system size. Extrapolation to 1/lniV = yields 
P max (oo) ~ 0.51. 



At b 1 the convergence to the fixed-point distribu- 
tion is much faster (Fig. |^b) . In this regime the condition 
for the scaling of the IPR distribution is N 1 / 2 jb 3> 1, as 
can be seen by comparing the relative magnitude of the 
IPR fluctuations at the critical point, [var^)] 1 / 2 ^^) — 
0.41/6 [see Eq. (pq)] with that in the Gaussian orthogonal 
ensemble, [var(P 2 )] 1/2 / \Pi) ^ IM/N 1 / 2 . Extrapolating 
the numerically found values of the relative variance to 
1/N — 0, we find good agreement with the theoretical 
prediction, as shown in Fig. [l0| Note a qualitative dif- 
ference in the approach to the fixed-point distribution 
at small and large b: while at b <C 1 the height of the 
distribution 7 3 (lnP2) decreases with N, approaching the 
limiting value from above, the behavior is opposite at 
fc> 1. 
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FIG. 10. Variance of P2 in the limit of large system sizes N 
for b = 100 (o), b = 8 (□) and b = 4 (O). The filled symbols 
denote the theoretical result (^) valid at 1 <C b <C V^iV, 
the solid line indicates the RMT limit, b 3> v~N, where 
[var(P 2 )] 1 / 2 /<P 2 ) « 1M/VN. 

Extrapolating the small-6 and large-6 results to the 
crossover range b ~ 1, we find simply the condition 
N 3> 1 for the proximity to the fixed-point distribution. 
Therefore, in the crossover regime b ~ 1 the limiting dis- 
tribution can be reached most easily. This is precisely 
what we observe in our numerical simulations. Figure 
||a, representing the evolution of the IPR distribution at 
6=1 demonstrates the almost perfect scale-invariance of 
the distribution with only very small deviations (less than 
3%) over a range of system sizes from N — 256 . . . 4096. 

Let us also comment on other types of "finite size ef- 
fects" that appear in numerical simulations. Numerically, 
it is impossible to perform an average at precisely a given 
value of the energy. Instead, one averages over an energy 
interval AE that one would like to choose as big as pos- 
sible in order to improve statistics. On the other hand, it 
is clear that AE should not be too big in order to avoid 
mixing of different critical theories. In our simulations 
we have chosen AE to be about 10% of the bandwidth. 
This value is still small enough, the corresponding vari- 
ation of the density of states v{E) being of the order of 
1%. 

Furthermore, the size 5 of the matrix ensemble the 
average was taken over is an important parameter in 
the simulations. Typical values we have used are: N — 
256,5' = 30000; N = 512,5* = 10000; N = 1024,5 = 
5000; N = 2048, 5 = 1000; N = 4096, 5 = 100. In some 
cases, like for the two- level correlation function i?2(s) or 
the full IPR distribution function, the ensemble has to 
be much larger. In these cases we give the numbers ex- 
plicitly in the figure caption. 



The compressibility \ can be expressed through the con- 
nected part i?2 ( s ) — ^2(5) — 1 of the critical level cor- 
relation function as follows: 



X 



dsR { 2 c) (s) . 



(63) 



Recently, it was arg ued in Ref. f| that Eq. (@) con- 
stitutes an exact relation between the spectral compress- 
ibility x an d the fractal dimension D2 ■ The derivation of 
(0) in is based on Dyson's idea of Brownian motion 
through the ensemble of Hamiltonians combined with 
some assumption of the decoupling of the energy-level 
and wave function correlations previously proposed in 
p9[ . While this decoupling has been proven to work up 
to three-loop order in the 1/g-expansion in 2D |39|], its 
applicability in the strong-coupling regime remained in 
the status of a conjecture. The PRBM model allows us 
to check the validity of the relation (Q) . Similarly to the 
IPR distribution function, the level correlation function 
can be calculated analytically in the two limits 6^-1 
and b <C 1 and numerically in the full range of b. 

In the b 3> 1 regime the two-level correlation function 
is obtained by an appropriate generalization of the ear- 
lier findings for the diffusive samples |^,^]; the results 
can be found in |2^^,|l^]. In particular, considering for 
simplicity the (3 = 2 ensemble at the band center, the 
level correlation function has the form 



sin 2 (7rs) [irs/Ab) 2 
(tts) 2 sinh 2 (7rs/46) 



(64) 



The correlation function ( |64| ) follows the RMT result 

^i\ s ) = ^( s ) — sin 2 (7rs)/(7rs) 2 up to the scale s ~ b 
(playing the role of the Thouless energy here), and then 
begins to decay exponentially. The spectral compress- 
ibility at b > 1 is given by [§0],[l5) 



III. LEVEL STATISTICS 

The two-level correlation function is defined in the 
usual way, 



1 



<"(£)> 



■(v(E + u/2)v(E-u/2)) , (61) 



where v(E) = N~ 1: Tr S(E — H) is the fluctuating density 
of states. At the critical point i?2(^) acquires a scale- 
invariant form (if considered as a function of s = u>/A, 
the frequency normalized to the mean level spacing A = 
l/N(v(E))) §^ff|. The distinct feature of the critical 
level statistics is a non-trivial value of the spectral com- 
pressibility < x < 1 characterizing the linear behavior 
of the variance of the number n{£ ) of levels in an energy 
window £ |30]]3|||3, 



var[n(£)] = x(n(£)) 



(n(£)) = I » 1 . (62) 



2tt/36 



6> 1 



(65) 



Comparing this with (g3|) , one finds [Q that the formula 
(0) is indeed satisfied to leading order in 1/6. 

We now turn to the opposite limit b <C 1. The evolu- 
tion equation for R%{u), r) can be written down in analogy 
with Eq. (ph: 



dR 2 (uj,r) _ 2_ 
' ~ N 



dr 



dr 



,IV-U[e- v2f2 l 2b2 2\V\ 
2?r b 



x [5(u> - 2V\J 1 + t 2 ) - 5(u) - 2Vt)] . (66) 



Equation ( p6| ) should be integrated over r from r = 
to A^/2 with the boundary condition i?2(w, 0) = 1; the 
result R2(u>, N/2) will then give the sought level corre- 
lation function. Evaluating the V^-integral in ( |66| ) and 
changing the variables to z = 2r/N and x = 2V/u, we 
get ats = w/ A>0 



11 



i r i 

s irz ■ 
dz — sin — 

7T0 



da: 



2 7 



x exp 



S X 2 7TZ 



(67) 



After some algebra, we find the level correlation function 
to be given by 



R { 2 c) (s) =<J(a)-crfc 



( 1*1 



(68) 



where erfc(x) = (2/y/iF) f°° exp(—t 2 )dt is the error func- 
tion and we have included the 8(s) contribution due to 
the self-correlation of the energy levels. Substitution of 
(E38|) in (63) yields the spectral compressibility 



X 



1-46 



&« 1 



(69) 



We see, therefore, that in the limit of small b the level 
repulsion is efficient in a narrow region |s| < b only, and 
the spectral compressibility tends to the Poisson value 
X = 1. Hence, formula (^), which would predict x — ► 1/2 
0, is violated. Similar violation of (0) is indi- 



at b 



catcd by numerical data for the tight-binding model in 
dimensions d > 4 [Q. Most likely it is never an exact 
relation, but rather an approximation valid in the weak- 
multifractality limit only. Strictly speaking, our results 
do not rule out the possibility that Eq. (0) is exact at b 
exceeding a certain value b c and breaks down at b < b c . 
However, we do not see a physical reason for such a qual- 
itative change being induced by the variation of b. 




FIG. 11. Two-level correlation function i?2(s) for two 
system sizes TV = 256 (o) and TV = 512 (□) at 
b = 0.1. The solid line indicates the theoretical re- 
sult (^). Number S of matrices in the ensemble: 
TV = 256, S = 3440512; TV = 512, 5" = 1418688. 

These results are fully supported by our numerical 
data. In particular, Fig. 11 represents the level corre- 



lation function ifefs) at b = 0.1 showing a nice agree- 
ment with Eq. $5%. Note that the finite size effects in 
the level correlation function at b <C 1 are much weaker 



than in the IPR distribution (Sec. IIP ). Indeed, the only 
assumption (besides b <C 1) used in the derivation of the 
evolution equation (^36j) is N 1, and then Eq. (p8|) is 
obtained without any further approximations. Therefore, 
in contrast to the IPR statistics, which reaches its fixed- 
point form at exponentially large N (the condition being 
folniV ^> 1), the level statistics acquires the IV- invariant 
form already at N > 1, see Figs, [ll], |l2[ At b > 1 the 
fixed-point condition is N b 2 (the same as for the IPR 
distribution, see Sec. HE). 
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FIG. 12. Variance of the number of levels {8n 2 } in a fixed 
energy interval as a function of the energy width of the in- 
terval parametrized by the mean level number (n) it con- 
tains. Traces correspond to b = 1 (open o: N — 4096, filled: 
TV = 2048), b = 0.25 (open □: TV = 4096, filled: TV = 2048) 
and b = 0.05 (open O: TV = 4096, O with dot: TV = 1024, 
filled: TV = 512). Statistical errors are typically of the order 
of the symbol size. The dashed line indicates the analytical 
prediction Eq. (|69|). 



To find numerically the spectral compressibility x, we 
plot the level number variance var[n(£ )] versus the aver- 
age (n(£ )) (Fig. |l^). The data show an extended plateau 
region in var[n(£ )]/{n(£ )}, determining \. The upper 
bound for this region is set by the matrix size TV, while 
the lower bound is ~ b (the value of the upper limit at 
which the integral ( |63|) saturates). We see that the data 
traces are independent of the system size TV (with ex- 
ception of the large- (n) cutoff determined by N) within 
the statistical errors. The numerically obtained spectral 
compressibility in the broad range of b is shown in Fig. [l^; 
in the large-& and small-& regions it agrees well with the 
corresponding analytical asymptotics. 
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FIG. 13. 



a function of 
> 1) to the 



Spectral compressibility \ as 
b: crossover from the "quasi-metallic" (6 
"quasi-insulating" (b <C 1) behavior. The lines indicate the 
analytical results for b S> 1 and i< 1, Eqs. j65| ) and (|69|). 

The above calculation is easily modified to the (3 = 2 
case, by replacing the measure of the ^-integration in 
( |66| ) according to (|60|). Performing the V- integral, we 
now get 



/f M, _ , /"' ^ s " „;..-> t: 



(«)=-! + / dz^ sin — 



o 



x exp 



2nb 2 



sm 



which yields the result 



R^ c \s) =S(s)-exp 



2tt& 2 



The spectral compressibility is thus equal to 
X ~ 1 - nV2b . 



(70) 



(71) 



(72) 



Again we see that the region of level repulsion shrinks 
in the limit 6^0, with the compressibility approaching 
the Poisson value x = 1. 

It is worth mentioning that the above results for the 
case [3 = 2 can also be obtained by exploiting the relation 
Jl2],[l5| between the PRBM model and a random matrix 
ensemble introduced by Moshe, Neuberger, and Shapiro 
p4j| . This mapping, which becomes exact in both limits 
b -C 1 and 6^1, relates the level statistics (but not the 
eigenfunction statistics!) of the two ensembles. On the 
other hand, the level correlation function of the ensem- 
ble of Moshe et al. can be calculated exactly in the case 
(3 = 2 and is in fact identical to the density correlation 
function of a ID non-interacting Fermi gas at a finite 
temperature Q]. Applying the results of Q], one ob- 
tains for the (3 = 2 PRBM ensemble precisely the results 
(H), © at b » 1 HQ and @, © at b < 1 @. 



IV. CONCLUSIONS 

In this paper, we have presented a detailed study of the 
statistics of eigenfunctions and energy levels in the family 
of the critical PRBM models. We have obtained analyt- 
ical results for the IPR distribution function, the mul- 
tifractal spectrum and the level correlation function in 
the two limits of weak and strong multifractality (6 > 1 
and b -C 1). The analytical results are fully supported 
by numerical simulations, which also have allowed us to 
explore the crossover region (b ~ 1). 

On the qualitative level, our main findings can be sum- 
marized as follows: 

1 . The distribution function of the IPR (normalized to 
its typical value P q yp ) is scale-invariant in the limit 
of large system size N. In other words, the dis- 
tribution function of the IPR logarithm, V(laP q ) 
shifts along the x-axis with increasing N, without 
changing its form and width. 

2. The scaling of P q yp with the system size defines 
the fractal exponent D q , which is a non-fluctuating 
quantity. 

3. The scale- invariant distribution V(z = P q /P q yp ) 
has a power-law tail oc z~ 1 ~ Xq . At sufficiently large 
q one finds x q < 1, and the average value (P q ) be- 
comes non-representative and scales with a differ- 
ent exponent D q ^ D q . 

4. The critical spectral statistics shows a crossover 
from a "quasi-metallic" (close-to-RMT) behavior at 
b ^S> 1 to a "quasi-insulating" (close-to-Poisson) one 
at b <C 1. In particular, the spectral compressibil- 
ity changes from to 1, thus violating the relation 
(0) in the strong-multifractality regime. 

Finally, it is worthwhile to comment on the extent of 
universality of the IPR distribution. Like the conduc- 
tance distribution or the level statistics |^l| , the IPR dis- 
tribution at criticality does depend on the system geom- 
etry (i.e., on the shape and on the boundary conditions). 
However, for a given geometry it is independent of the 
system size and of microscopic details of the model, and 
is an attribute of the relevant critical theory. 

Discussions with Y.V. Fyodorov, V.E. Kravtsov, 
L.S. Levitov, D.G. Polyakov, I. Varga and I.Kh. Zhareke- 
shev are gratefully acknowledged. This work was sup- 
ported by the SFB 195 der Deutschen Forschungsgemein- 
schaft. 
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